Perturbation method to model enamel caries progress 
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We develop a theoretical model of the carious lesion progress caused by acids diffusing into the 
tooth enamel from the dental plaque. The acids react with static hydroxyapatite, which leads to 
demineralization of the enamel, and consequently to the development of the carious lesion. The 
model utilizes the diffusion-reaction equations with one static and one mobile reactant where the 
reaction term is proportional to the product of concentrations of acids and of mineral. The changes 
of concentrations are calculated approximately by means of a perturbation method. The analytical 
approximate solutions are compared with the numerical ones and experimental data. 

PACS numbers: 87.10.+e, 82.39.Rt, 87.90.+y, 05.40.-a 



I. INTRODUCTION 



The carious lesion progress in dental enamel has been 
extensivel y st udied experimentall y [ ll, [2| 
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The experiments were per- 



formed in vitro with extracted teeth being demineralized 



in buffer solutions of or 



anic acids at different values of 
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These studies examined the time evolution of the depth 
of the carious lesion [E [E [IE [Hi 
mineral loss in the dental enamel 
rate of enamel dissolution 



20 



21 



_ 

;namel JlJ^JJJ, Jl7 , 



3C|, the 
1£|, the 
26j, the 



affect of various factors on the progress of caries such as 
concentrations and pH of buffer or inhibitors [IE HE E3 , 
the concentration profiles of hydroxyapatite (HA) (which 
is the main component of enamel) [I], IE H, H, HE HI, HE 
HI HE HE US [H, HE HE] . However, there are only a few 
attempts at theoretical description of the carious lesion 
progress. 

The formation of carious lesion of enamel starts when 
concentration of organic acids in the dental plaque 
reaches the sufficient value and pH of dental plaque low- 
ers below the appropriate point. Then, the organic acids 
diffuse in undissociated and/or dissociated form inward 
the enamel and react with the mineral to form soluble cal- 
cium and phosphate ions (or complexes) [E [IE HE H|[ . 
Thus, theoretical model should describe the diffusion of 
acids inside the enamel followed by the reaction with 
static hydroxyapatite. The diffusion-reaction equations 
for the system with one static and one mobile reactants 
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are usually chosen as 
dAjx.t) 

at 

dCjx.t) 

at 



D 



d 2 A(x,t) 
dx 2 

R(x,t), 



R(x,t), 



(1) 
(2) 



where A is the concentration of diffusing particles, D - 
their diffusion coefficient, C is the concentration of static 
reactant, and R denotes the reaction term, which is cho- 
sen within the mean field approximation as [32], [33| 



R(x,t) = kA m (x,t)C n {x,t), 



(3) 
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k is the reaction rate; the parameters n and m (which 
may be non-integer) are determined experimentally. The 
explicit solutions of Eqs. ([1]) and @ with the reac- 
tion term [Eg . ([3])] was found only for very few special 
cases [IE [34[- Thus, one is usually interested in find- 
ing some general characteristic functions of the diffusion- 
reaction system such as the time evolution of reaction 
front, the width of depletion zone etc. [35| . which are 
usually found by means of numerical methods of solv- 
ing the diffusion-reaction equations and computer simu- 
lations. To obtain approximate concentration profiles A 
and C simplifying methods are used, as for example the 
quasistatic [36| or perturbation [13, 0, HH ones. 

To theoretically describe the caries one often adopts 
assumptions which oversimplify the problem. Some au- 
thors consider only one of Eqs. (HJ and ^j. Additionally 
the reaction term R is taken in oversimplified form, where 
only one reactant is taken into account. For example, 
only Eq. |T]) was considered in [l7| with R = k(C s — C), 
where C s is a constant related to 'solubility' of the so- 
lute. In the reaction term was chosen in the form 
R(x,t) = k(x, t)A(x, t). The model based on Eq. ((U) 
only was studied in [l3| with R(x,t) = fcC n , and in [j] 
with R(x,t) = k(C s - C). Both Eqs. © and © were 
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considered in [U|4l|], where i?(x,i) = kA. We add that 
there were also considered the models based only on the 
diffusion equation without chemical reactions included 
through the 'effective' diffusion coefficient [rl El flU Et| . 
Very special diffusive theoretical models were used to de- 
scribe some caries characteristics such as temporal evo- 
lution of the amount of hydroxyapatite in the buffer 
solution released form the enamel during its dissolu- 
T3l 31 , l4ll. time evolution of the caries 
[TOnOTlMTll El, and concentration pro- 
files of fluoride @ and of hydrogen ions [4(| inside the 
enamel in the stationary state. As far as we know, there 
the theoretical concentrations of the mineral as a func- 
tion of time have not yet been obtained. 

Characteristic feature of carious lesion process is the 
creation of the surface layer in which the loss of hydrox- 
yapatite is relatively small. However, an appearance of 
this layer complicates theoretical models because the acid 
particles that diffuse through the layer do not chemically 
react with hydroxyapatie. The existence of a surface 
layer is not taken into account explicitly in the theoretical 
models mentioned above. 

In our paper we use both of Eqs. ([TJ) and ^ to describe 
the carious lesion process. The surface layer is also in- 
cluded in our considerations. We focus our attention on 
finding theoretical formulas of the concentration profiles 
of both reactants during the caries progress. To solve 
Eqs. ([H and (J2]) we use a perturbation method. The 
analytical approximate solutions are compared with the 
numerical ones of the diffusion-reaction equations and 
experimental data. 



II. CARIOUS LESION PROCESS 

The enamel consists of hydroxyapatite crystals 
Caio(P0 4 )(OH) 2 in 92 - 94% by weight H]. These 
crystals are organized in larger forms called prisms. The 
intercrystallinc and intcrprismatic spaces of enamel are 
filled with water EE HE]. Because of spaces between 
crystals and prisms, enamel is a microporous mate- 
rial @, i, |, EE EE EE, El- Besides HA the enamel 
includes inorganic factors, mostly fluoride and carbon- 
ate. In addition more then 40 trace elements can oc- 
cur in the tooth mineral. Organic matrix represents 
less then 1% by weight of the enamel [43|. The surface 
of dental enamel is covered by the dental plaque which 
mainly consists of microorganisms, saliva, leftovers and 
mucus. Oral microorganisms metabolize simple sugars 
coming from diet [E EE] to the organic acids (e.g., acetic 
or lactic). As we have mentioned previously, the forma- 
tion of carious lesion of enamel starts when concentra- 
tion of organic acids in dental plaque reaches sufficient 
value and pH of the dental plaque lowers below the ap- 
propriate point. Then the organic acids diffuse into the 
enamel @,H, H, El EE] • The acids can be transported in 
dissociated or undissociated form that depends on pH of 
the dental plaque @, i, EE El EE] . After achieving the 
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FIG. 1: The schematic view of the tooth enamel. The doted 
line represents the concentration of hydroxyapatite. 



enamel interior, acid reacts with the mineral according 
to the chemical formula [E EE, El El EE El 

Ca w {PO i ) 6 {OH) 2 + SH + -> 10C7a 2+ + 6H PO 2 ^ +2H 2 0, 



where the phosphate ions have an acidic form determined 
by the pH of the system. The products of reaction are 
inert for the caries progress Q . It is commonly accepted 
that the products of the reaction, calcium ions and phos- 
phate ions (or complexes) .are transported out of the 
enamel by the diffusion Q, i, 0, i, El EE El EE EH • 
There are two distinguishing stages in the formation of 
caries. In the first stage an apparently intact surface 
layer is created 0, i, EE El EE El, El (see Fig. [TJ. 
This is the layer where the loss of mineral is small in 
comparison to the content of the mineral in the sound 
enamel. Two possible mechanisms responsible for cre- 
ating the surface layer have been proposed. First, the 
appearance of protective agents prevents acids from dis- 
solving the enamel @, ©, EE El El HE E3- Second, 
the surface layer appears as a result of the combination 
of the dissolution and reprecipitation processes EE Ell • 
The thickness of this layer after reac hing a maximum 
value remains unchanged later on d EE]- Dissolution 
of the subsurface mineral (situated below the apparently 
intact surface layer) occurs in the second stage of the 
formation of carious lesion. The dissolution reaction 
takes place in a restricted region of the enamel called 
the reaction zone (the inner zone or the decalcification 
region). At the beginning the zone is placed right below 
the outer enamel surface but after exhausting the hydrox- 
yapatite in that reg ion the reaction zone moves inside 
the enamel [E IS LLJ] . Host factors involved in the caries 
process are as follows: composition and structure of the 
enamel [E E6l. p H and buffer concentration El E~~ 
kind of acid [id, EE, EJ , mineral content gradient 
composition of saliva, age of tooth, environment 
hygiene, etc. 
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III. MODEL 

The system under study is three-dimensional, but it 
is assumed to be homogeneous in the plane perpendicu- 
lar to the x axis, which is normal to the tooth surface. 
So, it can be treated as one-dimensional. The system 
described in the previous section is presented schemati- 
cally in Fig. [TJ This scheme is based on the qualitative 
descriptions of caries presented in the papers [H, |28| and 
on the plots of exp erimental concentration profiles given 

in 0, H 0, i, Ha El US US US [S3, Hi, US M- 

The border between the surface layer and the subsurface 
one is not sharp in realistic system, but for simplicity 
we assume that these regions are separated by the plane 
located at x = d. 

In the following we label the surface layer that does 
not change in time as / and the subsurface layer as II. 
Let us denote the concentration of hydrogen ions as A in 
the / region, B in the II region and the concentration 
of HA in the II region as C. We assume that the pure 
diffusion occurs in the region /, and there is the diffusion 
with chemical reactions in the region //. The ions dif- 
fuse in both regions with diffusion coefficient D. There 
arises a problem with a choice of the reaction term, as 
the parameters m and n occurring in Eq. ^ have not 
been unambiguously determined experimentally or the- 
oretically for the caries process. Beside we note that m 
and n depend on the chemical composition and structure 
of the tooth enamel, which is a personality trait, and 
on tooth environment, which seem to be uncontrollable 
in real systems. Here we assume that the reaction term 
is given by Eq. ([3]) with m = n = 1. This assumption 
simplifies the procedure of solving Eqs. {1} and ([2]). We 
presume that the case with m > 1 and/or n > 1 will 
not change significantly the qualitative form of the solu- 
tions (see for example [33| , where the cases of m = n = 1 
and m = n = 2 are considered). Thus, we choose the 
equations describing the caries as follows 



dA(x,t) _ D d 2 A(x,t) 
dt dx 2 ' 

dB(x,t) d 2 B{x,t) 



dt 
dC{x,t) 
Ft 



D- 



dx 2 

= -kB(x,t)C(x,t) 



(4) 

kB(x,t)C(x,t), (5) 
(6) 



At the initial time t = we assume that the enamel is 
free of acid and HA forms a homogeneous medium of the 
concentration Co- This assumption provides the initial 
conditions 



A(x,0) = B{x,0) = 0, 



C(x,0) = C . 



(7) 



(8) 



The dental plaque is the reservoir of the acid, so we as- 
sume that at the border of the enamel (x = 0) the con- 
centration of the acid Aq is constant. The thickness of 



the enamel is very large compared to the region where 
the acid concentration varies, so we treat the enamel as 
semi-infinite medium. It is rather obvious that the acid 
concentration and flux are continuous at the border be- 
tween the / and 77 regions. Thus we adopt the following 
boundary conditions 



A(0,t) = A ,A(d,t) = B(d,t), 



dA(x, t) 



dx 



dB(x,t) 



Ox 



B(po,t) = 0. 



(9) 



(10) 



(11) 



Because there is no exact analytic method to solve 
Eqs. ([4])- ©, we use the perturbation technique to find 
approximate solutions. 



IV. PERTURBATION METHOD 

To use the perturbation method we first transform 
Eqs. ([4l)- (fTTj) to the dimensionless form using the sub- 
stitutions: 



x = px s , t = rt B 



(12) 



where p and r denote the dimensionless position and 
time, x s and t s are constants of the dimension of space 
and time, respectively. Using Eq. (|12p we obtain the 
transport equations 



da(p,r) = d 2 a(p,T) 
dr dp 2 

db{p, t) _ d 2 b(p,r) 
dr dp 

dc(p, t) 



(13) 

, sb{p,r)c{p,T), (14) 



dr 

initial conditions 



erb(p,T)c(p,T), 



(15) 



o(p,0)=6(p,0)=0,c(p,0) = l, (16) 
and boundary conditions 

a(0, t) = l,a(p d , t) = b(p d , t), (17) 



da(p, t) _ db(p, t) 
Op lp=pd dp '" =pd ' 

b(oo, t) = 0, 



(18) 
(19) 



where 



a(p,r) 
b(p,r) 
c{p,r) 



A(px s ,Tt s ) 
A • 

BjfiXg, Tt s ) 

A 

C(pX s ,Tt s ) 

Co ' 



(20) 
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and 



A Q d 

G x s 



kC t s 



(21) 



(22) 



tions 



Equations (|13|) and (|T4|) are derived under condition that 
the parameters a; s and t s fulfill the relation 



'Dt a 



(23) 



Let us assume now that the parameter e is small 
enough (e <C 1)) to apply the perturbation method (with 
respect to this parameter) to find the concentration of 
HA as a function of space and time. 

Within the perturbation method the concentrations 
are given as 



n=0 
oo 

b(p,r) = J2 b n(p,r)s n , 

oo 



(24) 
(25) 
(26) 



Substituting Eqs. (|2"I ]) -([2l) ]l to Eqs. ([T5 ]) -(|Tg jl and compar- 
ing the functions of the same order with respect to the 
parameter e occurring at both sides of these equations, 
we get the equations of the zeroth order 



n=0 



dag(p,T) _ d 2 a {p, t) 

dr dp 2 
db (p,T) d 2 b {p,T) 



dr 
dr 

with the initial conditions 



dp 2 



0, 



ao(p, 0) = 0, 
6o(p,0) = 0, 
c (p,0) = 1, 



and boundary ones 



o (0,t) = 1, 
ao(Pd,r) = b (p d ,T), 
da (p,T) db (p,r) 



(27) 
(28) 
(29) 



(30) 
(31) 
(32) 



(33) 
(34) 

(35) 



dp lp=pd ~ dp lp=pd - 

For the order of n — 1, 2, 3, ... we get the following equa- 



da n {p,T) 

dr 
db n (p,r) 

dr 



dc u (p,r) 
dr 



d 2 a n (p, t) 

dp 2 
d 2 b n (p,r) 

dp 2 

n— 1 

-^2 b k(p, r)c n _ k -i(p,r) 



k=0 
n-1 

-r 

k=0 



^6 k (p, r)cn_ k _i(p,r) 



with the initial conditions 



O n (p,0) = 0, 

6 n (p,0) = 0, 

Cn(p,0) = 0, 



and the boundary conditions 



a n (0,r) = 0, 
a n (p d ,T) = b n (p d ,r) 



da n (p,r) 



dp 



db D (p,r) 



dp 



(36) 

(37) 
(38) 



(39) 
(40) 
(41) 



(42) 
(43) 

(44) 



We solve the above equations by means of the Laplace 
transform method [iij . Equations (|36p - (j38|) become 
more and more difficult to solve when the order of pertur- 
bation method grows. In practice, it is possible to obtain 
explicit solutions up to the first order. In the next sec- 
tion we find the exact solutions for ai, b\ and c\ where 
i = 0,1. 



V. ANALYTICAL AND NUMERICAL RESULTS 

Solving Eqs. ([2"7 P -([2l? l) with initial [Eqs. ([3Ti ]) -([S2 ]) ] and 
boundary [Eqs. (f33|) - (f35|) ] conditions we obtain 



a (p,T) = b (p,r) = crfc 



2^ 



c (p,T)=l, (45) 



where erfc(u) = exp(—r] 2 )dri is the complemen- 

tary error function. So, the zeroth order solutions take 
form of the pure diffusion without chemical reactions. 
Equations of first order are 



dai(p,r) _ d 2 ai{p, t) 

dr dp 2 
dh(p,T) 9 2 &i(p,r) 



dr 
dr 



dp 2 

-rb {p, t)c q {p,t) 



b (p,T)c (p,T), (46) 



When supplemented by the initial conditions ai(p, 0) 
b±(p, 0) = ci(p, 0) = and the boundary ones ai(0,r) 
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0,ax(p d ,r) = h(p d ,r) and ^f^| P=Pd - 



the solutions of Eqs. (|46|) are 

ai(p, r) = 



, f erfc ( ^JL 



-(2p d + p) A /-e 



-(2p d +p) 2 /4r 



{2 Pd -p)\j2 Pd - P 



2 J~~\ 2^ 



lU T+ (2 Pd + P )^ ei{ j2 Pd + P 



2 I V 2V? 

(2p d + p)J-e^+^^ 



+ (4 Pd -3p)J- e ^ 2 /4. 



T _ : 
IT 



(47) 



(48) 




FIG. 2: The concentrations of hydroxyapatite (increasing 
functions) and acid (decreasing functions) in the subsurface 
layer as a function of distance from the tooth surface for dif- 
ferent times given in the legend, all quantities are given in 
dimensionless units. The parameters are e = 1(T 5 and r = 1. 
The continuous lines with black symbols represent the ap- 
proximate solutions (|50|l - (|52p . the continuous lines with white 
symbols-the numerical solutions of Eqs. (|13[) -(115 [ ). 



( r+ £) erfc 


p 







(49) 

In Figs.(2Hl]we present the plots of the numerical solu- 
tions of Eqs. (P~5|) - (P~5]) and the dimensionless concentra- 
tions of acid 

a(p,T)=a (p,T)+eax(p,T), (50) 

for < p < p d , 

Hp, t ) = b o(p,r) +eb 1 (p,r), (51) 

for p > p d (the decreasing functions of p) and the di- 
mensionless concentration of HA (the growing functions 
ofp) 

C (P, T ) = co(p, t) + eci(p, t). (52) 

To use the numerical procedure we take the following 
approximations of the derivatives dc(p, t)/8t — [c(p, t) — 
c( P , t - At)}/ At, d 2 c(p, r)/dp 2 = [c(p + Ap, r) + c{p - 
Ap, t) — 2c(p, r)]/(Ap) 2 , calculations were performed for 
Ap = 0.25 and At = 0.01. The numerical solutions and 
approximate analytical ones are computed for few values 
of £ and r and for several values of time indicated in the 
legends of the plot; in all cases p d — 5. 

For e = 10~ 5 and r = 1 the numerical solutions (con- 
tinuous lines with white symbols) and approximate ana- 
lytical ones (continuous lines with black symbols) repre- 
senting the acid concentrations coincide with each other 




FIG. 3: The numerical and approximate analytical solutions 
as in Fig. [2]but for e = 10 -4 and r = 1. 



(Fig. [2|). Such a coincidence occurs for larger values of 
the parameters only for relatively small times (Figs. [3] 
and |4]) . Similar behavior is observed for the solutions 
representing the concentration of HA, but the difference 
between numerical and approximate analytical solutions 
are larger for the latter functions. Let us note that up 
to the first order approximation [Eqs. ([50]) - ([52)) ] the con- 
centrations of diffusive reactant depend explicitly on pa- 
rameter e whereas the concentrations of static reactant 
depend on the product of the parameters sr. When er is 
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FIG. 4: The description is the same as in Fig. [2] but for 
e = 10~ 4 and r = 10. 



FIG. 5: Theoretical concentration of HA (continuous lines) 
calculated for if = 5.5 h and At = 48 h and numerical solu- 
tions C(x,t) (dashed lines), separated points represent exper- 
imental data taken from [30]. 



relatively small, then the loss of mineral is noticeable for 
relatively large times. For example, the 5% loss of min- 
eral is observed after t w 6000 if er = 10~ 5 , after t « 1000 
if er = 10~ 4 , and after t w 100 if er = 10~ 3 . In the plots 
presented in Fig. [2HU we assume that er <C 1. Similar 
plots can be obtained for relatively large er for sufficiently 
small times. Of course, the above remarks concern the 
dimensionless quantities. To express the functions in the 
dimensional variables one needs to set the values of x s 
and t B obtained on the basis of experimental data. 



VI. COMPARISON WITH EXPERIMENTAL 
DATA 

The plots of HA shown in Fig. qualitatively agree 
with the experimental data. Fig. [2] is similar to Fig. 1 
from Fig. 1 from [1| and Figs. 3 A and 3 B from 0], 
Fig. [3] is like Fig. 3 from [yL Fig. 0] is in qualitative ac- 
cordance with Fig. 1 from [l[. Quantitative comparison 
of the model predictions with the experimental data ap- 
pears to be difficult. Experimental results obtained at 
very similar conditions are often very different from each 
other. The loss of HA achieves ~ 80% of the initial value 
within 14 days according to [26] but only ~ 40% of the 
initial value within 89 days as reported in the paper . 
This can be caused by individual features of each tooth. 
The amount of experimental data is too small to deter- 
mine the influence of the various factors. 

Despite the difficulties, we analyze the data taken from 
Fig. 1 presented in [3(| and from Fig. 2 presented in [29[ 
(C is given in special units g/cm 2 commonly used in scan- 
ning microradiography). The experimental concentration 
profiles extracted from these figures and theoretical func- 
tions are shown in Figs.[5]and[51 respectively. Taking into 
account how the experimental points are scattered, and 
including the accuracy of our readout, we estimated the 




x, mm 



FIG. 6: Experimental results and theoretical ones obtained 
for t{ — 17 h and At — 54 h; additional description is in the 
text. 



error to be 32 /mi for x and 8 mg/cm 2 for C(x,t) in Fig. [5] 
and 0.02 mm for x, 0.0096 g/cm 2 for C(x,t) in Fig. [6] 
These errors are marked in Figs. [5] and [HI The theoret- 
ical concentration profile expressed through dimensional 
variables is obtained from Eqs. (pT] ) -(f23 )) , (US]) , (T49]) . (j52j) . 
and it reads as 



C(x,t) = C Q -A C k 



x 2 \ f 
t + — ^ erfc 

2D J \2VDt 




-x 2 /4Dt 



(53) 



In our model we assume that the enamel is homoge- 
neous at the initial moment. This assumption is not al- 
ways fulfilled in real system. The initial concentrations 
are different for different tooth which is caused by indi- 
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vidual features of each tooth such as hygiene, environ- 
ment, etc. The experimental profiles corresponding to 
the first observation are also different from our initial 
condition ((5J) in the cases considered in our paper. Tak- 
ing into consideration the initial condition given by a 
specific function different from ([8|) leads to the complica- 
tions in the calculations and this can cause the pertur- 
bation method to be practically useless. Therefore, we 
assume that the initial condition can be chosen in the 
form given by Eq. ([5]), but we introduce a parameter tf 
which is interpreted as a time interval during which the 
demineralization process leads form homogeneous con- 
centration of HA to the 'initial concentration' presenting 
in the plots. Of course, the parameter t{ is not a 'real 
time' of demineralization process for a tooth situated in 
its natural environment; it represents a time of achiev- 
ing the initial concentration in the experimental system. 
Thus, the time of observation is given as 

t = t { + nAt, 

where n = 0, 1, 2, . . ., and At is the time step of successive 
concentration measurements. 

In Figs. [5] [5] the separated points represent the ex- 
perimental data, the continuous lines are the approxi- 
mate solutions given by Eq. (|53[) and the dashed lines 
are the numerical solutions of Eqs. ([4])-([6]). The fitting 
parameters which ensure the best matching of approxi- 
mate solutions (|53|) to experimental data are if, D and 
the product A$k. The most natural choice of the pa- 
rameter t s is t s = l/(Aok) what leads to er = 1 [see 
Eqs. (Ell) and @] and x s = y/D/(A k). In the papers 
[29l [301 ] the initial concentration of A is given in the units 
of mol/cm 3 and there is rather impossible to express A 
in the units used in the plots. Since the parameters A 
and k are not known, to obtain the numerical solutions 
we use the equations written in the dimensionless form 
(fT3|) -(fl"5" |) with er = 1. The values of e and r are not 
determined unambiguously by experimental data. The 
numerical calculations performed for different values of 
the parameters show that the functions c(p, r) obtained 
for any values of e under condition that e < 10 -3 and 
r = 1/e do not differ significantly from each other. The 
reason is that the time evolution of function c(p, r) de- 
pends explicitly on the product er [Eq. (fT5")) ] and its de- 
pendence on e is only manifested by the function b(p, t), 
which for small e is very close to the function bo {p, r) . 
So, we chose the parameters e = 1/r = 10~ 3 for nu- 
merical calculations. In considered cases the subsurface 
layer is rather thin in both of them. The height of this 
layer reaches about 30 per cent of the initial values of the 
concentration of HA at the second measurement and it 
decreases within time. So, we conclude that its influence 
on the concentration profiles (especially for long times) is 
small. For considered cases we assume that the width of 
the layer is very small comparing to the subsurface layer 
and we take pd = 2 to calculations. After finding the 
solutions of dimensionless equations, we transform them 
to the dimension form according to the formula (|20p . 



In Fig. [5] the functions were calculated for the param- 
eters Co = 87 mg/cm 2 and time interval At — 48 h 
taken from [30]. The fit parameters are found as A^k = 
1.5 x 10" 4 1/s, D = 2.2 x lO" 10 cm 2 /s and t { = 5.5 h. In 
Fig. [6] we obtained the functions for Co = 0.116 g/cm 2 , 
D = 9.9 x 10~ 13 cm 2 /s, A k = 1.5 x 10~ 3 1/s and 
t[ = 17 h. As far as we know the exact values of the pa- 
rameters k and D are unknown for the substances under 
considerations. Moreover, it is clear that these values de- 
pend on individual features of a tooth. Thus, we can not 
conclusively compare the obtained value of D with the 
experimental one mentioned in the literature. However, 
we note that the parameter takes 'typical' values. For ex- 
ample it was reported in 0] that D oc 10 _8 ~10 -13 cm 2 /s 
for several substances diffusing in the tooth enamel. 

We conclude this section by saying that the model 
qualitatively describes the data. The quantitative com- 
parison is not fully conclusive but the fit is quite satis- 
factory. 



VII. FINAL REMARKS 

We propose a theoretical model of the carious lesion 
progress, which is based on physically well-motivated 
diffusion-reaction equations. In contrast to the previ- 
ous theoretical models of caries, the reaction term de- 
pends on the concentrations of diffusing acid and of 
enamel mineral. To approximately solve the diffusion- 
reaction equations we adopt a perturbation method. 
The concentration profiles presented in our paper are 
in qualitative agreement with the experimental data 

el s i n m m m mm, he m m, mm except m 

the region where the concentration of hydroxyapatite is 
the smallest. In the theoretical model the concentration 
of hydroxyapatite goes to zero in the region located just 
beyond the surface layer. In a real system (schematically 
shown in Fig. [I} a non-zero concentration of enamel is 
always observed. This is caused by occurrence of various 
impurities, such as fluoride, in the enamel. These impuri- 
ties prevent a complete loss of the enamel. However, tak- 
ing the impurities into considerations significantly com- 
plicates the model. We are mostly interested in modeling 
of caries in the region, where changes of HA concentra- 
tion are the largest because therein the caries progress is 
most intensive and the limit of caries is located in this 
region. Thus, we do not take into account the enamel 
impurities in the considered model. Analyzing the ex- 
perimental concentration profiles of HA presenting in the 
paper cited above, one can see that the surface layer is 
often created in the form differing from the one assumed 
in our model, for example the HA concentration decrease 
in time in this layer. Moreover, the initial concentra- 
tion of HA is not given by Eq. ([5]). In our opinion the 
model considered in our paper can be used to describe 
the caries in such situations. We observed that the re- 
duction of the HA concentration in surface layer does not 
noticeably change the solutions of the diffusion-reaction 
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equations. Introduction of the parameter t a allows one 
to find the solution of the equations corresponding to the 
initial concentration in the experiment. 

A quantitative comparison requires returning to the di- 
mensional variables x and t. This is possible if one knows 
values of the parameters occurring in Eqs. ((I])-©. As far 
as we know, the values of these parameters are unknown 
for teeth studied in the papers presenting the measured 
concentrations. Qualitative comparison is also difficult 
because of individual features of each tooth caused by 
various factors such as age, environment, hygiene, etc., 
which result in the different quantitative characteristics 
of each tooth. Thus, repeatability of the experimental 
results is weak. The changes of concentration of the acid 
are controlled mainly by the parameter e only, whereas 
the loss of hydroxyapatite mostly depends on the prod- 
uct of the parameters er. We note that the parameter 
r = Aq/Cq is controlled in experiments in vitro because 
Aq is the acid concentration of the buffer solution. 

The perturbation method was already used to solve 



the diffusion-reaction equations for the homogeneous sys- 
tem consistin g o f two mobile initially separated reac- 
tants [37l. l38l. l39l| . where, similarly to our paper, up to 
the first order solutions were obtained. However, the first 
order solutions in [371 . l38l . |39I | were calculated only within 
long time approximation. Let us note that in our case of 
one static and one mobile reactant, we found the exact 
solutions up to the first order within the perturbation 
method. 
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